/*
 * Copyright 1993-2006 NVIDIA Corporation.  All rights reserved.
 *
 * NOTICE TO USER:
 *
 * This source code is subject to NVIDIA ownership rights under U.S. and
 * international Copyright laws.
 *
 * This software and the information contained herein is PROPRIETARY and
 * CONFIDENTIAL to NVIDIA and is being provided under the terms and
 * conditions of a Non-Disclosure Agreement.  Any reproduction or
 * disclosure to any third party without the express written consent of
 * NVIDIA is prohibited.
 *
 * NVIDIA MAKES NO REPRESENTATION ABOUT THE SUITABILITY OF THIS SOURCE
 * CODE FOR ANY PURPOSE.  IT IS PROVIDED "AS IS" WITHOUT EXPRESS OR
 * IMPLIED WARRANTY OF ANY KIND.  NVIDIA DISCLAIMS ALL WARRANTIES WITH
 * REGARD TO THIS SOURCE CODE, INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY, NONINFRINGEMENT, AND FITNESS FOR A PARTICULAR PURPOSE.
 * IN NO EVENT SHALL NVIDIA BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL,
 * OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
 * OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE
 * OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE
 * OR PERFORMANCE OF THIS SOURCE CODE.
 *
 * U.S. Government End Users.  This source code is a "commercial item" as
 * that term is defined at 48 C.F.R. 2.101 (OCT 1995), consisting  of
 * "commercial computer software" and "commercial computer software
 * documentation" as such terms are used in 48 C.F.R. 12.212 (SEPT 1995)
 * and is provided to the U.S. Government only as a commercial end item.
 * Consistent with 48 C.F.R.12.212 and 48 C.F.R. 227.7202-1 through
 * 227.7202-4 (JUNE 1995), all U.S. Government End Users acquire the
 * source code with only those rights set forth herein.
 */

/* Matrix multiplication: C = A * B.
 * Host code.
 */

// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

// includes, project
#include <cutil.h>

// includes, kernels
#include "MN_cuda_kernel.cu"
#include "MN_cuda.h"


////////////////////////////////////////////////////////////////////////////////
// declarations, forward

extern "C"
void computeGold(float*, const float*, const float*, unsigned int, unsigned int, unsigned int);

void CompactMatrix(int n, int p, float *S_h, 	float *Sf_ptr_rm, float *Sf_idx_rm, float *Sf_data_rm, 
												float *Sr_ptr_rm, float *Sr_idx_rm, float *Sr_data_rm, 
												float *Sf_ptr_cm, float *Sf_idx_cm, float *Sf_data_cm, 
												float *Sr_ptr_cm, float *Sr_idx_cm, float *Sr_data_cm);

void ReadFile(float* M, char* file_name, unsigned int size);
void mnSeq(int n, int p, float *vf, float *vr, float *vt, float *a, float *dm, float *m, float *kf, float *kr, float *S, float epsilon, float tolerance, float t);

void ComputeOnDevice(float *S_h, float *kf_h, float *kr_h, float *m_h, float *vt_h, int n, int p);//, float epsilon, float tolerance, float t);

/*float S_test[48] = {
				1,-1,-1,0,0,0,0,0,
				0,0,2,-1,-1,0,0,0,
				0,1,0,0,-1,0,1,0,
				0,1,0,0,0,0,-1,-1,
				0,0,0,0,1,-1,0,0,
				0,0,0,1,0,-1,0,0
			};
			
float kf_test[8] = {1,1,1,1,1,1,1,1};
float kr_test[8] = {0,1,1,1,1,0,1,0};
float m_test[6] = {1,1,1,1,1,1};
*/
//float epsilon_test = 0.01;
//float tolerance_test = 0.00001;
//float t_test = 0.001;

float epsilon = 0.01;
float tolerance = 0.00001;
float t = 0.001;

//int n_test = 8;
//int p_test = 6;

////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////
int main(int argc, char** argv) {

	
	//srand(52);
	
	float*  S;
	float*  Kf;
	float*  Kr;
	float*  m;
	float*  vfc;
	float*  vrc;
	float*  vtc;
	float*  ac;
	float*  dmc;
	float*  vt;

	unsigned int n,p;
	uint GPUTime,CPUTime;

	cutCreateTimer(&CPUTime);
    	cutCreateTimer(&GPUTime);
    	cutResetTimer(CPUTime);
    	cutResetTimer(GPUTime);	
	
	srand(2012);
	
	if(argc == 6 ) 
	{
		// Allocate and read in matrices from disk
		unsigned int k;
		int* params = NULL; 
		unsigned int data_read = 2;
		cutReadFilei(argv[1], &params, &data_read, true);
		
		n = params[0];
		p = params[1];
		k=n*p;

		printf("n = %d, p = %d\n", n, p);

		S  = (float*) malloc(k*sizeof(float));
		Kf = (float*) malloc(n*sizeof(float));	
		Kr = (float*) malloc(n*sizeof(float));
		m  = (float*) malloc(p*sizeof(float));
		vt = (float*) malloc(n*sizeof(float));
		vfc = (float*) malloc(n*sizeof(float));
		vrc = (float*) malloc(n*sizeof(float));
		vtc = (float*) malloc(n*sizeof(float));
		ac = (float*) malloc(p*sizeof(float));
		dmc = (float*) malloc(n*sizeof(float));
		
		cutFree(params);
		
		(void)ReadFile(S, argv[2], k);
		(void)ReadFile(Kf, argv[3], n);
		(void)ReadFile(Kr, argv[4], n);
		(void)ReadFile(m, argv[5], p);
		
		for(int i = 0; i < 42; i++)
			printf("S[%d] = %3.4f\n", i, S[i]);
		
		//ComputeOnDevice(S.elements, kf.elements, kr.elements, m.elements, vt.elements, S.width, S.height, params[2], params[3], params[4]);
		
		//FreeMatrix(&S);
    		//FreeMatrix(&kf);
    		//FreeMatrix(&kr);
		//FreeMatrix(&m);
	}
	else //if we can't read the input, then run test
	{
			//float *vt;
			//vt = (float*) malloc(n_test*sizeof(float));
	       printf("unable reading parameter file\nExecuting test\n");
		   //ComputeOnDevice(S_test, kf_test, kr_test, m_test, vt, n_test, p_test, epsilon_test, tolerance_test, t_test);
			
	}

	cutStartTimer(CPUTime);
	mnSeq(n,p,vfc,vrc,vtc,ac,dmc,m,Kf,Kr,S,epsilon,tolerance,t);
	cutStopTimer(CPUTime);

	cutStartTimer(GPUTime);
    	ComputeOnDevice(S, Kf, Kr, m, vt, n, p);
	cutStopTimer(GPUTime);


	for(int i = 0; i < n; i++)
	{
	    printf("vtc[%d] = %f, vt[%d] = %f\n", i, vtc[i], i, vt[i]);
	}
	
	printf("Time for CPU: %f ms\n",cutGetTimerValue(CPUTime));
	printf("Time for GPU: %f ms\n",cutGetTimerValue(GPUTime));    

	return 0;
}

////////////////////////////////////////////////////////////////////////////////
//! Run a simple test for CUDA
////////////////////////////////////////////////////////////////////////////////
void ComputeOnDevice(float *S_h, float *kf_h, float *kr_h, float *m_h, float *vt_h, int n, int p)//, float epsilon, float tolerance, float t)
{

	vt_h = (float*) malloc(n*sizeof(float));
	
	float *S_d;
	float *kf_d;
	float *kr_d;
	float *vf_d;
	float *vr_d;
	float *vt_d;
	float *dm_d;
	float *m_d;
	float *a_d;
	
	
    	cudaMalloc((void**)&S_d, n*p*sizeof(float));
	cudaMalloc((void**)&kf_d, n*sizeof(float));
	cudaMalloc((void**)&kr_d, n*sizeof(float));
	
	cudaMalloc((void**)&vf_d, n*sizeof(float));
	cudaMalloc((void**)&vr_d, n*sizeof(float));
	cudaMalloc((void**)&vt_d, n*sizeof(float));
	cudaMalloc((void**)&dm_d, (p+2)*sizeof(float));
	cudaMalloc((void**)&m_d, p*sizeof(float));
	cudaMalloc((void**)&a_d, p*sizeof(float));
	
	
	cudaMemcpy(S_d, S_h, n*p*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(kf_d, kf_h, n*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(kr_d, kr_h, n*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(m_d, m_h, p*sizeof(float), cudaMemcpyHostToDevice);
	
	cudaMemset(a_d,1,p*sizeof(float));
	cudaMemset(dm_d,1,(p+2)*sizeof(float));
	cudaMemset(vt_d,1,n*sizeof(float));
	cudaMemset(vf_d,1,n*sizeof(float));
	cudaMemset(vr_d,1,n*sizeof(float));
	
	
	
	// Setup the execution configuration
	dim3 dim_grid, dim_block;
	cudaError_t cuda_ret;
	
	// Launch the device computation threads!
	/* Kernel Call */
	/*cudaDeviceProp dev_prop; 
	cudaGetDeviceProperties( &dev_prop, 0);
	printf("Max Threads per block = %d\n", dev_prop.maxThreadsPerBlock);
	printf("Shared Memory per block = %d\n", dev_prop.sharedMemPerBlock);
	printf("Max threads per multiprocessor = %d\n", dev_prop.maxThreadsPerMultiProcessor);
	printf("MultiProcessor count = %d\n", dev_prop.multiProcessorCount);*/
	
	//int TILE_WIDTH = 16;
	
	dim_block.x = (n>p) ? n : p;
	dim_block.y = 1; 
	dim_block.z = 1;
	dim_grid.x = 1;
	//if(MATRIX_WIDTH % dim_block.x != 0) dim_grid.x++;
	dim_grid.y = 1;
	//if(MATRIX_HEIGHT % dim_block.y != 0) dim_grid.y++;
	dim_grid.z = 1;

	int max_iterations = 1000; //this limits how many iterations to run
	
	unsigned int computeOnDevice;
	cutCreateTimer(&computeOnDevice);

	/////////////
	//unoptimized kernel
	/////////////
	cutStartTimer(computeOnDevice);
	
	MN_cuda_Kernel<<<dim_grid,dim_block>>>(max_iterations, n, p, vf_d, vr_d, vt_d, a_d, dm_d, m_d, kf_d, kr_d, S_d, epsilon, tolerance, t);
	cuda_ret = cudaDeviceSynchronize();
    	if(cuda_ret != cudaSuccess) printf("Unable to launch/execute kernel");
		
	cutStopTimer(computeOnDevice);
	printf("Total time on kernel 1 %f\n", cutGetTimerValue(computeOnDevice));
	
	/////////////
	//kernel2 with register, shared memory, and math optmization
	////////////
	int shared_mem_bytes = 2*n+2*p+2; //2*n for vf and vr, p for m, p+2 for dm with the 2 for reduction with p odd
	cutResetTimer(computeOnDevice);
	cutStartTimer(computeOnDevice);
	
	MN_cuda_Kernel2<<<dim_grid,dim_block,(shared_mem_bytes)*sizeof(float)>>>(max_iterations, n, p, vt_d, m_d, kf_d, kr_d, S_d, epsilon, tolerance, t);
	cuda_ret = cudaDeviceSynchronize();
    	if(cuda_ret != cudaSuccess) printf("Unable to launch/execute kernel");
		
	cutStopTimer(computeOnDevice);
	printf("Total time on kernel 2 %f\n", cutGetTimerValue(computeOnDevice));
	
	///////////
	//kernel3 with matrix compaction
	//////////
	float Sf_ptr_rm[p+1]; //Sf array of pointers to beginning of row
	float Sf_idx_rm[n*p]; //Sf column indices in row major
	float Sf_data_rm[n*p]; //Sf non zero data elements in row major
	float Sr_ptr_rm[p+1]; //Sr array of pointers to beginning of row
	float Sr_idx_rm[n*p]; //Sr column indices in row major
	float Sr_data_rm[n*p]; //Sr non zero data elements in row major
	float Sf_ptr_cm[n+1]; //Sf array of pointers to beginning of column
	float Sf_idx_cm[n*p]; //Sf column indices in column major
	float Sf_data_cm[n*p]; //Sf non zero data elements in column major
	float Sr_ptr_cm[n+1]; //Sr array of pointers to beginning of column
	float Sr_idx_cm[n*p]; //Sr column indices in column major
	float Sr_data_cm[n*p]; //Sr non zero data elements in column major
	CompactMatrix(n, p, S_h, 	Sf_ptr_rm, Sf_idx_rm, Sf_data_rm, 
									Sr_ptr_rm, Sr_idx_rm, Sr_data_rm, 
									Sf_ptr_cm, Sf_idx_cm, Sf_data_cm, 
									Sr_ptr_cm, Sr_idx_cm, Sr_data_cm);
									
	float *Sf_ptr_rm_d;
	float *Sf_idx_rm_d;
	float *Sf_data_rm_d;
	float *Sr_ptr_rm_d;
	float *Sr_idx_rm_d;
	float *Sr_data_rm_d;
	float *Sf_ptr_cm_d;
	float *Sf_idx_cm_d;
	float *Sf_data_cm_d;
	float *Sr_ptr_cm_d;
	float *Sr_idx_cm_d;
	float *Sr_data_cm_d;
	
	cudaMalloc((void**)&Sf_ptr_rm_d, (p+1)*sizeof(float));
	cudaMalloc((void**)&Sf_idx_rm_d, (n*p)*sizeof(float));
	cudaMalloc((void**)&Sf_data_rm_d, (n*p)*sizeof(float));
	cudaMalloc((void**)&Sr_ptr_rm_d, (p+1)*sizeof(float));
	cudaMalloc((void**)&Sr_idx_rm_d, (n*p)*sizeof(float));
	cudaMalloc((void**)&Sr_data_rm_d, (n*p)*sizeof(float));
	cudaMalloc((void**)&Sf_ptr_cm_d, (n+1)*sizeof(float));
	cudaMalloc((void**)&Sf_idx_cm_d, (n*p)*sizeof(float));
	cudaMalloc((void**)&Sf_data_cm_d, (n*p)*sizeof(float));
	cudaMalloc((void**)&Sr_ptr_cm_d, (n+1)*sizeof(float));
	cudaMalloc((void**)&Sr_idx_cm_d, (n*p)*sizeof(float));
	cudaMalloc((void**)&Sr_data_cm_d, (n*p)*sizeof(float));
	
	cudaMemcpy(Sf_ptr_rm_d, Sf_ptr_rm, (p+1)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sf_idx_rm_d, Sf_idx_rm, (n*p)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sf_data_rm_d, Sf_data_rm, (n*p)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sr_ptr_rm_d, Sr_ptr_rm, (p+1)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sr_idx_rm_d, Sr_idx_rm, (n*p)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sr_data_rm_d, Sr_data_rm, (n*p)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sf_ptr_cm_d, Sf_ptr_cm, (n+1)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sf_idx_cm_d, Sf_idx_cm, (n*p)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sf_data_cm_d, Sf_data_cm, (n*p)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sr_ptr_cm_d, Sr_ptr_cm, (n+1)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sr_idx_cm_d, Sr_idx_cm, (n*p)*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(Sr_data_cm_d, Sr_data_cm, (n*p)*sizeof(float), cudaMemcpyHostToDevice);
	
	shared_mem_bytes = 2*n+2*p+2; //2*n for vf and vr, p for m, p+2 for dm with the 2 for reduction with p odd
	cutResetTimer(computeOnDevice);
	cutStartTimer(computeOnDevice);
	
	MN_cuda_Kernel3<<<dim_grid,dim_block,(shared_mem_bytes)*sizeof(float)>>>(max_iterations, n, p, 
																	vt_d, m_d, kf_d, kr_d, 
																	Sf_ptr_rm_d, Sf_idx_rm_d, Sf_data_rm_d, 
																	Sr_ptr_rm_d, Sr_idx_rm_d, Sr_data_rm_d, 
																	Sf_ptr_cm_d, Sf_idx_cm_d, Sf_data_cm_d, 
																	Sr_ptr_cm_d, Sr_idx_cm_d, Sr_data_cm_d, 
																	epsilon, tolerance, t);
	cuda_ret = cudaDeviceSynchronize();
    	if(cuda_ret != cudaSuccess) printf("Unable to launch/execute kernel");
		
	cutStopTimer(computeOnDevice);
	printf("Total time on kernel 3 %f\n", cutGetTimerValue(computeOnDevice));
		
		
	// Read P from the device
	cudaMemcpy(vt_h, vt_d, n*sizeof(float), cudaMemcpyDeviceToHost);
	
	for(int i = 0; i<n; i++) {
			printf("vt[%d] = %f\n",i,vt_h[i]);
		}
		
	cudaFree(S_d);
	cudaFree(kf_d);
	cudaFree(kr_d);
	cudaFree(vf_d);
	cudaFree(vr_d);
	cudaFree(vt_d);
	cudaFree(dm_d);
	cudaFree(m_d);
	cudaFree(a_d);
	//free(vt_h);
	
	
}

//Compact matrix
void CompactMatrix(int n, int p, float *S_h, 	float *Sf_ptr_rm, float *Sf_idx_rm, float *Sf_data_rm, 
												float *Sr_ptr_rm, float *Sr_idx_rm, float *Sr_data_rm, 
												float *Sf_ptr_cm, float *Sf_idx_cm, float *Sf_data_cm, 
												float *Sr_ptr_cm, float *Sr_idx_cm, float *Sr_data_cm)
{
	//calculate Sf row major
	int nElements = 0;
	int nElementsPrev = 0;
	int nElementsRow = 0;
	memset(Sf_ptr_rm, 0, (p+1)*sizeof(float));
	memset(Sf_idx_rm, 0, (n*p)*sizeof(float));
	memset(Sf_data_rm, 0, (n*p)*sizeof(float));
	
	for(int i = 0; i<p; i++) {
		for(int j = 0; j<n; j++) {
			if(S_h[n*i+j]<0) {
				Sf_data_rm[nElements] = -1*S_h[n*i+j];
				Sf_idx_rm[nElements] = j;
				nElements++;
			}
			Sf_ptr_rm[i] = nElementsRow;
		}
		nElementsRow = nElementsRow + nElements - nElementsPrev;
		nElementsPrev = nElements;
	}
	Sf_ptr_rm[p] = nElements;
	
	//calculate Sr row major
	nElements = 0;
	nElementsPrev = 0;
	nElementsRow = 0;
	memset(Sr_ptr_rm, 0, (p+1)*sizeof(float));
	memset(Sr_idx_rm, 0, (n*p)*sizeof(float));
	memset(Sr_data_rm, 0, (n*p)*sizeof(float));
	
	for(int i = 0; i<p; i++) {
		for(int j = 0; j<n; j++) {
			if(S_h[n*i+j]>0) {
				Sr_data_rm[nElements] = S_h[n*i+j];
				Sr_idx_rm[nElements] = j;
				nElements++;
			}
			Sr_ptr_rm[i] = nElementsRow;
		}
		nElementsRow = nElementsRow + nElements - nElementsPrev;
		nElementsPrev = nElements;
	}
	Sr_ptr_rm[p] = nElements;
	
	//calculate Sf column major
	nElements = 0;
	nElementsPrev = 0;
	nElementsRow = 0;
	memset(Sf_ptr_cm, 0, (n+1)*sizeof(float));
	memset(Sf_idx_cm, 0, (n*p)*sizeof(float));
	memset(Sf_data_cm, 0, (n*p)*sizeof(float));
	
	for(int i = 0; i<n; i++) {
		for(int j = 0; j<p; j++) {
			if(S_h[n*j+i]<0) {
				Sf_data_cm[nElements] = -1*S_h[n*j+i];
				Sf_idx_cm[nElements] = j;
				nElements++;
			}
			Sf_ptr_cm[i] = nElementsRow;
		}
		nElementsRow = nElementsRow + nElements - nElementsPrev;
		nElementsPrev = nElements;
	}
	Sf_ptr_cm[n] = nElements;
	
	//calculate Sr column major
	nElements = 0;
	nElementsPrev = 0;
	nElementsRow = 0;
	memset(Sr_ptr_cm, 0, (n+1)*sizeof(float));
	memset(Sr_idx_cm, 0, (n*p)*sizeof(float));
	memset(Sr_data_cm, 0, (n*p)*sizeof(float));
	
	for(int i = 0; i<n; i++) {
		for(int j = 0; j<p; j++) {
			if(S_h[n*j+i]>0) {
				Sr_data_cm[nElements] = S_h[n*j+i];
				Sr_idx_cm[nElements] = j;
				nElements++;
			}
			Sr_ptr_cm[i] = nElementsRow;
		}
		nElementsRow = nElementsRow + nElements - nElementsPrev;
		nElementsPrev = nElements;
	}
	Sr_ptr_cm[n] = nElements;
}									

void ReadFile(float* M, char* file_name, unsigned int size)
{
	cutReadFilef(file_name, &M, &size, true);
	
}

void mnSeq(int n, int p, float *vf, float *vr, float *vt, float *a, float *dm, float *m, float *kf, float *kr, float *S, float epsilon, float tolerance, float t)
{

int i, j;
float fwdStoic[n * p];  //substrate stoichiometric values in forward direction - Sf
float revStoic[n * p];  //substrate stoichiometric values in reverse direction - Sr
unsigned int timeStep;

//forward and reverse stoichiometric matrices
for(i = 0; i < n * p; i++)
{
       if(S[i] < 0)
       {
           fwdStoic[i] = -1 * S[i];
           revStoic[i] = 0.0f;
       }
       else if(S[i] > 0)
       {
           fwdStoic[i] = 0.0f;
           revStoic[i] = S[i];  
       }
       else
       {
           fwdStoic[i] = 0.0f;
           revStoic[i] = 0.0f;  
       }
}

for(timeStep = 0; timeStep < 1000; timeStep++)
{
// loop through the fluxes and calculate the fwd and reverse flux vectors
    for(i = 0; i < n; i++)
    {
        vf[i] = kf[i];
        vr[i] = kr[i];
        for(j = 0; j < p; j++)
        {
            vf[i] = vf[i] * powf(m[j],fwdStoic[(j * n) + i]); // indices are reversed in matlab code
            vr[i] = vr[i] * powf(m[j],revStoic[(j * n) + i]);
        }
    }

    // calculate the metabolic concentration change

    for(i = 0; i < p; i++)
    {
        dm[i] = 0;
        for(j = 0; j < n; j++)
        {
                dm[i] += S[(i * n) + j] * (vf[j] - vr[j]);
        }
    }

    // Compute new metabolite concentrations
    for(i = 0; i < p; i++)
    {
        m[i] = m[i] + (a[i] * dm[i] * t);
    }

    for(i = 0; i < p; i++)
    {
        // if concentration is negative, set it to a very small value.
        if(m[i] < 0)
        {
            m[i] = 0.00001;
        }
    }

    // Compute new volume scale factor
    for(i = 0; i < p; i++)
    {
        a[i] = epsilon + ( 1 / (abs(dm[i]) + epsilon));
    }
}



// Calculate total flux
for(i = 0; i < n; i++)
{
   vt[i] = vf[i] - vr[i];
}

}
